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ABSTRACT 

We present numerical simulations of axisymmetric, magnetically driven relativistic 
jets. Our special-relativistic, ideal-MHD numerical scheme is specifically designed to 
optimize accuracy and resolution and to minimize numerical dissipation. In addition, 
we implement a grid-extension method that reduces the computation time by up to 
three orders of magnitude and makes it possible to follow the flow up to six decades 
in spatial scale. To eliminate the dissipative effects induced by a free boundary with 
an ambient medium we assume that the flow is confined by a rigid wall of a prescribed 
shape, which we take to be z oc r a (in cylindrical coordinates, with a ranging from 1 to 
3). We also prescribe, through the rotation profile at the inlet boundary, the injected 
poloidal current distribution: we explore cases where the return current flows either 
within the volume of the jet or on the outer boundary. The outflows are initially cold, 
sub-Alfvenic and Poynting flux-dominated, with a total-to-rest-mass energy flux ra- 
tio fi ~ 15. We find that in all cases they converge to a steady state characterized 
by a spatially extended acceleration region. The acceleration process is very efficient: 
on the outermost scale of the simulation as much as ~ 77% of the Poynting flux has 
been converted into kinetic energy flux, and the terminal Lorentz factor approaches its 
maximum possible value (Too ~ /u). We also find a high collimation efficiency: all our 
simulated jets (including the limiting case of an unconfined flow) develop a cylindrical 
core. We argue that this could be the rule for current-carrying outflows that start 
with a low initial Lorentz factor (To ~ I). Our conclusions on the high acceleration 
and collimation efficiencies are not sensitive to the particular shape of the confining 
boundary or to the details of the injected current distribution, and they are quali- 
tatively consistent with the semi-analytic self-similar solutions derived by Vlahakis 
& Konigl. We apply our results to the interpretation of relativistic jets in AGNs: wc 
argue that they naturally account for the spatially extended accelerations inferred in 
these sources (r^ > 10 attained on radial scales R > 10 17 cm) and are consistent with 
the transition to the matter-dominated regime occurring already at R > 10 16 cm. 
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1 INTRODUCTION 

There is strong evidence for relativistic motions in jets that 
emanate from active galactic nuclei (AGNs). In particular, 
apparent superluminal speeds /3 app (in units of the speed of 
light c) as high as ~ 40 have been measured for radio compo- 
nents on (projected) scales of ~ 1 — 10 pc in the blazar class 
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of sources (e.g. Ijorstad et al]|200ll) . Ijorstad et ail (|2005l ) 
used a method based on a comparison between the time 
scale of flux density decline and the light-travel time across 
the imaged emission region to relate /3 app to the bulk Lorentz 
factor r of the outflow; they inferred that the Lorentz fac- 
tors of blazar jets lie in the range ~ 5 — 40, with the ma- 
jority of quasar components having T ~ 16 — 18 and with 
BL Lac objects po ssessing a more uniform V distribution. 
ICohen et alj l|2007l ) recently reached a similar conclusion on 
the basis of probability arguments, inferring that roughly 
half the sources in a flux density-limited, beamed sample 
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have a value of V close to the measured /3 apP - They further 
deduced that the maximum Lorentz factor in their sample 
of 119 AGN jets is ~ 32 , close to the v alue of ~ 4 inferred 
for the jets observed bv lJorstad et al.1 (|200ll . \2Q0$) . 

The presence of relativistic bulk motions in blazar 
jets has been independently indicated by measurements 
of rapid variations in the total and polarized fluxes (e.g . 
Hartman et~aDl200ll : iRebillot et ai1l2006l : iBach et alJl200d 
Villata et al.ll2006f ). There is also evidence that the relativis- 
tic speeds persist to large scales. For example, apparent su- 
perluminal component motions have been measured in the 
3C 120 jet ou t to projected dista nces from the source of at 
least 150 pc (j Walker et al.ll200ll ). and it has been argued 
that the spectral properties of the heads of extended (up 
to several hundred kiloparsecs) jets can be explained in the 
context of a relativistic flow that is decelerated to subrela- 
tivistic speeds at the termination shock that advances into 
the ambient medium (jGeorganopoulos fc Kazanasll2003h . 

The main source of power of AGN jets is the rota- 
tional energy of the central supermassiv e black hole (e.g. 
lLovelacdfl976l : Blan dford fc Znaieklll977l) and/or its accre- 
tion disk (e.g. iBlandfordl 119761 1. The naturally occurring 
low mass density and hence high magnetization of black- 
hole magnetospheres suggests that the relativistic jets orig- 
inate directly from the black-hole ergosphere, whereas the 
disk surface launches a slower, possibly nonrelativistic wind 
that surrounds and confines the highly relativistic flow. This 
picture finds support in recent numerical simulatio ns (e.g 
iMcKinnev fc Gammiel I200I ; |Pe Villiers et al.ll2005h . How- 
ever, this issue is far from settled and one cannot rule out 
the inner part of accretion disk as a base of a relativistic out- 
flow fe.g. IVlahakis fc Ko nigl 2004). The theory of relativis- 
tic, magnetically driven jets from black holes (and neutron 
stars) predicts highly magnetized flows, with the Poynting 
flux dominating the total energy output. At the jet emis- 
sion site this energy has to be transferred to particles. This 
transfer may have the form of magnetic di ssipation (e.g. 
IBlandfordl |2002| ; iLvutikov fc Blandford! I2003T I but the still 
commonly held view is that the Poynting energy is first con- 
verted into bulk kinetic energy and only subsequently chan- 
neled into radiation through shock s and other dissipativ e 
waves (e.g. iBlandford fc Reeslll974l ; iBegelman et al.lll984h . 
The jet radiative efficiency is still one of the key debatable 
issues in the theory of Poynting-dominated outflows. In prin- 
ciple, slow magnetic dissipation in an expanding jet may also 
facilitate the g radual conversion of Poynt ing flux into bulk 
kinetic energy (|Drenkhahn fc Spruitll2002l ). 

When the inertia of the plasma is negligibly small its 
dynamics is well described by the approximat ion of force- 
free e l ectrodynamics (or magnetodynamics; e. g. iKomissarovl 
120021 ; iKomissarov. Barkov fc Lvutikovl l2007h . The equa- 
tions of magnetodynamics are much simpler than those of 
magnetohydrodynamics (MHD), which is what prompted 
their application to the study of the magnetic accel- 
eration of relativistic jets (e.g. I Blandfordl 1 19761 . |2002| ; 
iNaravan. McKinnev fc Farmerll2007l ). The solutions of these 
equations describe trans-Alfvenic flows whose drift veloc- 
ity approaches the speed of light at infinity, the location of 
the fast critical point in these models. However, within this 
framework it is impossible to account for the conversion of 
Poynting flux into plasma kinetic energy and to study the 
issue of the acceleration efficiency efficiency. 



The next simplest approximation that can be used to 
address the issue of Poynting-to-kinetic energy conversion 
is ideal MHD. In this case one can obtain exact semi- 
analytic solutions, although, because of the complexity of 
the problem, this can only be done when the system pos- 
sesses a high degree of symmetry. T his approach was pio- 
neered bv lBlandford fc Pavnd (| 19821 ), who constructed non- 
relativistic semi-analytic solutions for steady-state, cold, 
self-similar (in the spherical radial coordinate) disk out- 
flows. These solutions w ere generalized to the (special) rela- 
t iv istic MHD regime bv lLi. Chiueh fc Begelmanl (|l992l ) and 
by IContopoulosM 19941). T hey were further investigated by 
IVlahakis fc Konigll (|2003^ lbh. who also considered the ef- 
fect of thermal forces during the early phases of the accel- 
eration!]] Solutions with s imilar properties were derived in 
iBeskin fc Nokhrinal (|2006l ) by linearizing about a force-free 
solution for a paraboloidal field geometry. 

A key property of the relativistic solutions derived in 
the aforementioned studies is the extended nature of the 
acceleration region: the bulk of the (poloidal) acceleration 
is effected by magnetic pressure gradients (associated with 
the azimuthal magnetic field component) and takes place be- 
yond the classical fast-magn etosonic point (a singular point 
of the Bernoulli equation). IVlahakis fc Konigll (|2003al ) in- 
terpreted this b ehaviour (whic h was d ubbed the "magne tic 
nozzle" effect bv lLi et aljfl99i see also [Camenzin d 1989) m 
terms of the distinction between the classical and the modi- 
fied fast-magnetosonic surfaces fe.g. lBogovalov|[r997h . They 
pointed out that the latter surface, which is the locus of the 
fast-magnetosonic singular points of the combined Bernoulli 
and trans-field (or Grad-Shafranov) equations, is the true 
causality surface (or "event horizon") for the propagation 
of fast waves when the shape of the field lines is obtained 
from the solution of the trans-field equation (with the classi- 
cal surface playing this role only when the shape of the flux 
surfaces is predetermined). They argued that, in this case, 
the acceleration continues all the way to (and possibly even 
past) the modified fast-magnetosonic surface, which can lie 
well beyond the classical one0 Another general property 
of the cold MHD solutions is that, in the current-carrying 
regime (where the poloidal components of the current den- 
sity and the magnetic field are antiparallel) they collimate 
(asymptotically) to cylinders. Furthermore, the asymptotic 
Lorentz factor corresponds to a rough equipartition between 
the Poynting and kinetic energy fluxes. 

The continuation of the acceleration process beyond 
the classical fast-magnetosonic surface is evidently a gen- 
eral characteristic of steady-state M HD solutions that ap - 
plies also to nonrelativistic jets (e.g. IVlahakis et al-l feoOO). 
This behaviour should, however, be more clearly discerned 
in observations of relativistic flows, where the proper speed 

1 IVlahakis fc Konigll <2003ah focused on flows whose initial 
poloidal velocity component is sub-Alfvenic, corresponding to the 
poloidal magnetic field component dominating the azimuthal field 
compon ent at the top of the disk, whereas IVlahakis fc Konigll 
(2003b) discussed the super- Alfvenic case in which the azimuthal 
component is dominant at the base of the flow. In this paper we 
only consider outflows of the first type. 

2 In the radia ll y self -similar solutions presented in 
IVlahakis fc Konigll (2003a), the modified fast-magnetosonic 
surface formally lies at an infinite distance from the origin. 
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T/3 can increase by a large factor between the classical and 
the modified singular surfaces. In contrast, the magnetic ac- 
celeration of non-relativistic flows is almost complete at the 
classical fast point. This striking difference has a very sim- 
ple origin. For nonrelativistic flows the criticality condition 
at the classical fast-magnetosonic point implies equiparti- 
tion between the magnetic energy and the kinetic energy of 
poloidal motion. The kinetic energy can therefore increase 
by at most a factor of 2 beyond this point. However, rel- 
ativistic flows remain magnetically dominated at the fast- 
magnetosonic point, which means that there is an ample 
remaining supply of magnetic energy that can be used for 
flow acceleration downstream of this point (e.g. iKomissarovl 
|2004) . 

In the case of AGNs there have indeed been indica- 
tions from a growing body of data that the associated rela- 
tivistic jets undergo the bulk of their acceleration on scales 
that are of the order of those probed by very-long-baseline 
radio interferometry. In one line of reasoning, the absence 
of bulk-Comptonization spectral signatures in blazars has 
been used to infer that jet Lorentz factors > 10 are only 
attained on scales > 10 cm (Sikoraetal ] I2005T ). There 
have also been explicit inferences of component accelera- 
tion based on radio proper motion and X-ray e mission mea- 
surem ents for the jets in the quasa rs 3C 345 (|Unwin et alj 
1 19971 ) and 3C 279 jpiner et al.ll2003T ). Extended acceleration 
in the 3C 345 jet has been independently indicated by the 
higher apparent speed s of jet components loca ted further 
away from the nucleus l|Lobanov fc Rola nd 2005) and by the 
observed luminosity varia tions of the moving components 
ilLobanov fc Zensusj|l999l ). Similar effects in other blazars 



(e.g. iHoman et al.ll200lf ) suggest that parsec-scale accelera- 
tion to relativistic speed s may be a common feature of AGN 
jets. IVlahakis fc Konigll l|2004h argued that these observa- 
tions are most naturally interpreted in terms of magnetic 
driving and employed self-similar relativistic jet solutions to 
generate model fits to the 3C 345 data in support of this 
conclusion. 

While the semi-analytic solutions have been useful in 
indicating basic properties of the magnetic acceleration pro- 
cess and in providing valuable clues to the interpretation of 
the observational data, more general solutions are needed 
to confirm these results and to gain a fuller understanding 
of the generation of relativistic jets in AGNs. In particu- 
lar, numerical simulations are needed to find out whether 
the self-similar model captures the essential properties of 
outflows that obey realistic boundary conditions and that 
are not required to be in a steady state. Among the ques- 
tions that such simulations could answer are: (1) Do disk 
outflows in fact approach a steady state, and, if they do, 
is that state stable? (2) Is the acceleration indeed generally 
extended, and to what extent does the asymptotic state of 
the self-similar solutions approximate the far-field behaviour 
of more realistic outflows? (3) Do any new traits emerge 
when the restrictions imposed by the self-similarity assump- 
tion are removed? Of particular interest is the question of 
the ability of the magnetic driving mechanism to accelerate 
outflows to high Lorentz factors with high efficiency over 
astrophysically relevant distance scales. Another important 
question is whether highly relativistic flows can be strongly 
collimated by purely magnetic stresses. There have been lin- 
gering doubts over these issues in the literature (see Sec- 



tion [57T|, and although they have already received tentative 
answers, a full numerical study could help to settle them 
once and for all. 

Although there have already been several reported sim- 
ulations of the formation of jets in black-hole accretion flows 
using relativistic (in fact, generaZ-relativistic) MHD codes, 
so far they have provided only partial answers to the above 
questions. The existing calculations indicate that magnetic 
acceleration indeed operates over several decades in radius 
and can accelerate jets to relativistic speeds. However, the 
extended nature of the acceleration typically results in the 
bulk Lorentz factor reaching only a small fraction of its po- 
tential asymptotic value by the time the simulation is ter- 
minated. For example, in the longest jet simulated to date, 
which extended to ~ 10 4 times the gravi tational radius r g 
of the central black hole (|McKinnevll2006l ). the Lorentz fac- 
tor on the largest computed scale was ~ 10, which is just 
~ 10 -2 — 10~ 3 of the estimated asymptotic value. This im- 
pressive simulation deals with an extremely complex system 
of which the jet is only one component, the other being 
the black hole, the accretion disk, the disk corona, the low- 
speed "wall jet," and their surroundings. Ultimately, this is 
the kind of simulation one wants to carry out in order to 
fully understand the dynamics of AGN outflows. However, 
they are also very challenging from the computational point 
of view. One major concern ins that, in view of the extended 
nature of its acceleration, the jet is particularly vulnerable to 
numerical diffusion and dissipation. These numerical effects 
may partly explain why the quantity Too in the above-cited 
paper, which is the same as our n (equation I15|l and should 
be a field-line constant, in fact decreases by about one order 
of magnitude along a mid-leve l field line in the simulated jet 
(see Fig. 7 in lMcKinnevll200"6h . 

In this paper we address the above questions through 
numerical simulations specifically designed for investigat- 
ing the key aspects of ideal-MHD acceleration of relativis- 
tic jets. In the first place, we u se a numerical sc heme 
based on a linear Riemann solver l|Komissarovl Il999l ^ that 
does not need a large artificial diffusion for numerical 
stability. This distinguishes it from most other schemes 
for relativistic MHD, including those that are based on 



HLL , KT and similar flux prescriptions (e.g . Gammi e et al 



2006 



2003] : iDuez et all 120051: iKoide et al. I Il999l; lAnninos et al 



Shibata fc Sekuguchil 120051 : lAnderson et al 

200' 



2006a 



bel Zanna et al. I 120031 ; lAnton et al.l 120061 ). Simple one- 
dimensional tests suggest that this should lead to a no- 
ticeably greater accuracy in two-dimensional problems that 
involve stat ionary flows that are aligned with the computa- 
tional grid (Komissarov 2006). Secondly, instead of studying 
jet propagation through some ambient medium, we consider 
the case of a flow in a funnel with solid walls. This allows us 
to avoid the errors that would otherwise be caused by nu- 
merical mass diffusion and dissipation at the interface. Fi- 
nally, we employ elliptical (or spherical) coordinates adapted 
to the chosen paraboloidal (or conical) shape of the funnel. 
This allows us to have the jet well resolved everywhere (us- 
ing a fixed number of grid points across the funnel) and to 
benefit from the close alignment of the flow with the compu- 
tational grid. These careful measures in conjunction with a 
grid-extension method enable us, for the first time, to track 
the acceleration and collimation processes to their comple- 
tion. 
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We describe the basic equations in Section [2] and the 
numerical calculations in Section [3] The simulation results 
are presented in Section [4] and discussed in Section [5] We 
summarize in Section [6] 



2 BASIC EQUATIONS 

Since most of the acceleration takes place far away from 
the source, we assume that the space-time is flat. Moreover, 
the flow is described in an inertial frame at rest relative to 
the source. In this case we can write the system of ideal 
relativistic MHD as follows. The continuity equation 



-gpu 



) + a i ( v /= s/ow i ) = o J 



(i) 



where p is the rest mass density of matter, u v is its 4- velocity, 
and g is the determinant of the metric tensor; the energy- 
momentum equations 



(2) 



where T KU is the total stress-energy-momentum tensor; the 
induction equation 



(l/c)d t {B l )+e l]k d J {E k ) = Q, 



(3) 



where e,js = y/j^ijk is the Levi-Civita tensor of the absolute 
space (ei23 = 1 for right-handed systems and 6123 = — 1 for 
left-handed ones) and 7 is the determinant of the spatial part 
of the metric tensor (yy = gtj); the solenoidal condition 



0. 



(4) 



The total stress-energy- momentum tensor, T KV , is a 
sum of the stress-energy momentum tensor of matter 



T ( m ) =wu u /c +pg 



(•>) 



where p is the thermodynamic pressure and w is the enthalpy 
per unit volume, and the stress-energy momentum tensor of 
the electromagnetic field 



F 



*F» a -^(F^F aP )g K 



(6) 



where F VK is the Maxwell tensor of the electromagnetic field. 
The electric and magnetic field are defined as measured by 
an observer stationary relative to the spatial grid, which 
gives 

1 

2* 

and 



B* = ^e ijk F jk 



(7) 



E, = F it . (8) 
In the limit of ideal MHD 

Ei = -e m v d B k /c, (9) 

where v l — u 1 /u is the usual 3- velocity of the plasma. 
We use an isentropic equation of state 

p = Qp s , (io) 

where Q =const and s = 4/3. Since we are interested in the 
magnetic acceleration of cold flows, we make Q very small, 
so the gas pressure is never a dynamical factor. This rela- 
tion enables us to exclude the energy equation from the in- 
tegrated system. However, the momentum equation remains 



intact, including the nonlinear advection term. Therefore, if 
the conditions for shock formation were to arise, our calcu- 
lation would capture that shock0 



2.1 Field-line constants 

The poloidal magnetic field is fully described by the az- 
imuthal component of the vector potential, 

i,8Aa, 



B 1 = 



ij<t> ' 



(11) 



For axisymmetric solutions — ^ /2ir, where ^(x 1 ), 
the so-called magnetic flux function, is the total magnetic 
flux enclosed by the circle x 1 =const (x % being the coor- 
dinates of the meridional plane). Stationary and axisym- 
metric ideal MHD flows have 5 quantities that propagate 
unchanged along the magnetic field lines and thus are func- 
tions of ^ alone. These are k, the rest-mass energy flux per 
unit magnetic flux; il, the angular velocity of magnetic field 
lines; I, the the total angular momentum flux per unit rest- 
mass energy flux; /i, the total energy flux per unit rest-mass 
energy flux; and Q, the entropy per particle. For cold flows 
(Q — 0, w = pc 2 ) we have 



k = 



I = 



pu p 
B v 



VpB* 
r B v 



2ixkc 



+ ru 



and 



T(l + a) 



(12) 

(13) 
(14) 

(15) 



where u p = Tv p is the magnitude of the poloidal compo- 
nent of the 4-velocity, B v is the magnitude of the poloidal 
component of the magnetic field, r is the cylindrical radius, 

/ = JrB* (16) 

is the total electric current flowing through a loop of radius 
r, a is the ratio of the Poynting flux to the matter (kinetic 
plus rest-mass) energy flux, and 

r ^-2^ ^ 

is the Poynting flux per unit rest-mass energy flux. (Here 
and in the rest of the paper we use a hat symbol over vector 
indices to indicate their components in a normalized coordi- 
nate basis.) From equation (|15|l it follows that the Lorentz 
factor r cannot exceed p,. 



3 NUMERICAL SIMULATIONS 

To maintain a firm control over the jet's confinement and 
to prevent complications related to numerical diffusion of 

3 Since entropy is fixed the compression of our shocks would be 
the same as for continuous compression waves. This gives a higher 
jump in density for the same jump in pressure than in a proper 
dissipative shock. Fortunately, we do not need to contend with 
this issue in practice as shocks do not form in our simulations. 
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the dense nonrelativistic plasma from the jet's surroundings, 
we study outflows that propagate inside a solid funnel of 
a prescribed shapeQ Specifically, we consider axisymmetric 
paraboloidal funnels 



where z and r are the cylindrical coordinates of the funnel 
wall. This suggests the utilization of a system of coordinates 
in which the funnel wall is a coordinate surface. For a conical 
jet (a = 1) we use spherical coordinates, whereas for jets 
with a > 1 we employ elliptical coordinates {£, n, <j}}, where 

-l/a 



and 



V z 

a 



(18) 



(19) 



(see Appendix lAl for details) 

We use a Godu nov-type numerica l code based on the 
scheme described in iKomissarovl Jl999). To reduce numeri- 
cal diffusion we applied parabolic reconstruction instead of 
the linear one of the original code. Our procedure, in brief, 
was to calculate minmod-averaged first and second deriva- 
tives and use the first three terms of the Taylor expansion for 
spatial reconstruction. This simple procedure has resulted 
in a noticeable improvement in the solution accuracy even 
though the new scheme is still not 3rd-order accurate be- 
cause of the non-uniformity of the grid. 

The grid is uniform in the £ direction (the polar an- 
gle direction when we use spherical coordinates), where in 
most runs it has a total of 60 cells. To check the conver- 
gence, some runs were repeated with a doubled resolution. 
The cells are elongated in the r\ direction (the radial direc- 
tion when we use spherical coordinates) , reflecting the elon- 
gation of the funnel. For very elongated cells we observed 
numerical instability, so we imposed an upper limit of 40 on 
the length/width ratio. 

To speed up the simulations, we implemented a 
sectioning of the computation al grid as described in 
iKomissarov fc Lvubarskvl (|2004l ). In each section, which is 
shaped as a ring, the numerical solution is advanced using a 
time step based on the local Courant condition. It is twice as 
large as the time step of the adjacent inner ring and twice 
as small as the time step of the adjacent outer ring. This 
approach is particularly effective for conical flows but less 
so for highly collimated, almost cylindrical configurations. 

4 In real astrophysical systems, the shape of the boundary is de- 
termined by the spatial distribution of t he pressure or the density 
of the confining ambient medi um (e.g. iBlandford fc Reed 11974 : 
lKonig]|[T983; lKomissarovlll994) . The effective ambient pressure 
distributions implied by the adopted funnel shapes are consid- 
ered in Section l5.2l 

5 The equations are dimensionalized in the following manner. The 
unit of length, L, is such that r\i = 1 L, where the subscript i refers 
to the inlet boundary. The unit of time is T = L/c. The unit of 
mass is M = L 3 Bq/4ttc 2 , where Bo is the dimensional magnitude 
of the rj component of magnetic field at the inlet (so the dimen- 
sionless magnitude of at the inlet is V 4ir). In applications, 
L is the length scale of the launch region (e.g. the radius of the 
event horizon if the jet originates in a black hole), T is the light 
crossing time of that region and Bo is the typical strength of the 
poloidal magnetic field at the origin. 



3.1 Boundary conditions 

3.1.1 Inlet boundary 

We treat the inlet boundary, rji = 1, as a surface of a per- 
fectly conducting rotator and consider two rotation laws, 



n = fi 

and 

n = n [i-3(£/£j 



- 2(£/& 



(20) 



(21) 



where £j marks the jet boundary. The angular velocity pro- 
file is directly related to the distribution of the return elec- 
tric current in the jet (see equation 1281 below). In fact, the 
current is driven by the electric field associated with the ro- 
tating poloidal field, and charge conservation requires the 
circuit to eventually close. In the case of constant Q the re- 
turn current flows over the jet boundary. For the rotation 
law 1)210 it is distributed over the jet body as a volume cur- 
rent, the current changing sign at £ ~ £j/2. Thus, we cover 
the two generic types of electric current distribution. The 
solid-body rotation law provides a very good description of 
the behaviour of magnetic field lines that thread the horizon 
of a black hole. This choice is therefore entirely appropri- 
ate for the black-hole theory of relativistic AGN jets. The 
differential rotation law is more suitable to the accretion- 
disk theory, although it admittedly does not correspond to 
a realistic velocity field (which is hard to model given the 
limitations of our numerical technique). 

The condition of perfect conductivity allows us to fix the 
azimuthal component of the electric field and the r\ compo- 
nent of the magnetic field: 



(22) 



(23) 



(24) 



= 0, — Bo at 77 = ra . 
From the first of these conditions we derive 

and (using equation 1 13 p 



B 7 ! 

We have also experimented with non-uniform distributions 
of the magnetic field, in particular with B v decreasing with 
£. The results were not significantly different as the field 
distribution downstream of the inlet underwent a rapid re- 
arrangement that restored the transverse force balance. 

On the assumption of a cold (i.e. zero thermal energy) 
jet, the flow at the inlet boundary is necessarily super-slow- 
magnetosonic. This means that both the density and the 
radial component of the velocity can be prescribed some 
fixed values: 

p = po , v v = v P0 . 

In the simulations we used v pa = 0.5 c, which was a choice 
of convenience. On the one hand, this value is sufficiently 
small to insure that the flow at rj = 1 is sub-Alfvenic and 
hence that the Alfven and fast-magnetosonic critical sur- 
faces are located downstream of the inlet boundary. On the 
other hand, it is large enough to promote a rapid settle- 
ment to a steady state (keeping in mind that the speed of 
a steady-state flow remains constant along the symmetry 
axis). Because of the sub-Alfvenic nature of the inlet flow, 
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Figure 1. Model CI. Left panels show log 10 Yp (colour), where Yp is the jet density as measured in the laboratory frame, and magnetic 
field lines. Right panels show the Lorentz factor (colour) and the current lines. The thick solid line in the top-left panel denotes the 
surface where the flow becomes supcrfast in the r\ direction. The top panels show the solution for the first grid sector, whereas the bottom 
panels show the combined solution for the second and third grid sectors. 



we cannot fix the other components of the magnetic field 
and the velocity - they are to be found as part of the global 
solution. Following the standard approach we extrapolate 
and B^ from the domain into the inlet boundary cells. 
We then compute and from equations (|23|l and (|24|) . 

In the case of differential rotation the magnitude of the 
angular velocity is chosen in such a way that the Alfven 
surface of the jet is near the jet origin, its closest point being 
located at a distance of ~ 1.5 times the initial jet radius 
from the inlet surface. In the case of solid-body rotation 
the Alfven surface almost coincides with the light cylinder, 



whose radius r; c = c/fi is only 50% larger than the initial 
jet radius. 

The inlet density is chosen so that all jets have very 
similar values of /n and a. In particular, for the models with 
uniform f2 we have fj, maK ~ 18, and for the models with 
non-uniform SI we have ^t ma x — 12. 



3.1.2 Other boundaries 

The computational domain is always chosen to be long 
enough for the jet to be super-fast-magnetosonic when it 
approaches the outlet boundary r\ = r\ . This justifies the 
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Figure 2. Same as Fig.Q] but for Model C2. 



use of radiative boundary conditions at this boundary (i.e. 
we determine the state variables of the boundary cells via 
extrapolation of the domain solution) . 

At the polar axis, £ = 0, we impose symmetry boundary 
conditions for the dependent variables that are expected to 
pass through zero there, 

/(-0 = -/«)• 

These variables include B^, B^, and . For other vari- 
ables we impose a "zero second derivative" condition, 

d 2 f/df = o, 

which means that we use linear interpolation to calculate 
the values of these variables in the boundary cells. 

We do this in order to improve the numerical represen- 



tation of a narrow core that develops in all cases as a result 
of the magnetic hoop stress. Within this core the gradients 
in the £ direction are very large and the usual zero-gradient 
condition, /(— £) = /(f), results in increased numerical dif- 
fusion in this region. We have checked that this has a no- 
ticeable effect only on the axial region and that the global 
solution does not depend on which of these two conditions 
is used. 

At the wall boundary, f = f j, we use a reflection condi- 
tion, 

/(& + AO = -/(& - AO , 

for B^ and ir and a zero-gradient condition for all other 
variables. 
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Figure 3. Same as Fig. [T] but for Model A2. 



3.2 Initial setup 

The initial configuration corresponds to a non-rotating, 
purely poloidal magnetic field with approximately constant 
magnetic pressure across the funnel. The plasma density 
within the funnel is set to a small value so that the out- 
flow generated at the inlet boundary can easily sweep it 
away. In order to speed this process up the 77 component of 
velocity inside the funnel is set equal to 0.7 c, whereas the £ 
component is set equal to zero. 



3.3 Grid extensions 

The inner rings of the grid, where the grid cells are small 
and so is the time step, are the computationally most inten- 



sive regions of the simulation domain. If we kept computing 
these inner rings during the whole run then we would not be 
able to advance very far from the jet origin. Fortunately, the 
trans-sonic nature of the jet flow allows us to cease compu- 
tations in the inner region once the solution there settles to 
a steady state. To be more precise, we cut the funnel along 
the ^-coordinate surfaces into overlapping sectors with the 
intention of computing only within one sector at any given 
time, starting with the sector closest to the inlet boundary. 
Once the solution in the "active" sector settles to a steady 
state we switch to the subsequent sector, located further 
away from the inlet. During the switch the solution in the 
outermost cells of the active sector is copied into the cor- 
responding inner boundary cells of the subsequent sector. 
During the computation within the latter sector these in- 
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Figure 4. Flow constants fc(*&), O('I') and fJ.(^) at different distances from the source for models CI (top row), C2 (middle row) and A2 
(bottom row). In all cases the solid lines show the constants at the injection surface (t; = 1). The deviations from these values further 
downstream are due to a gradual accumulation of numerical errors. For models CI and C2 the dashed line corresponds to r) = 10 2 , the 
dash-dotted line to r} = 10 3 , and the dotted line to r] = 10 4 . For model A2 the dashed line corresponds to r\ = 2 X 10 2 , the dash-dotted 
line to rj = 2 X 10 3 , and the dotted line to r\ = 2 X 10 4 . 



ner boundary cells are not updated. Surely, this procedure 
is justified only when the flow in a given sector cannot com- 
municate with the flow in the preceding sector through hy- 
perbolic waves, and thus we need to ensure that the Mach 
cone of the fast-magnetosonic waves points outward at the 
sector interfaces. This condition can be written as 



where c 2 = spc 2 /w (see equation 1 10[) and b = (B 2 — E 2 ) 1 ^ 2 
is the magnetic field magnitude in the fluid frame (see Ap- 
pendix [B] for details). In the cold limit this reduces to 



(V/ c ) 2 -( 



v /c 



1 1/2 



Cf . 



(26) 



(rv) ; 



+- 



\ 4:TTW I C 



+ 



c 2 /c 2 Attw/c 2 



- cl/c 2 



+ 



>o, 



(25) 



where c/ = 6/(47rp + b /c ) ' is the isotropic fast speed in 
the fluid frame. Thus, the jet has to be super-fast in the 77 
direction at the sector interfaces. In Figs.[TH3]the location of 
the surface where v v = Cf is shown by a thick solid line: to 
the right of this line v v > Cf. One can see that the transition 
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Figure 5. Development of the line current in model C2. 
The figure shows the azimuthal magnetic field at r\ = 
200, 400, 800, 1600, 3200 and 6400 (from bottom to top). The 
solution at r] = 6400 is plotted as squares. 



to the super-fast regime occurs well inside the first sector. 
(Note that when v v > c/ the inequality 1261 is satisfied.) 

In these simulations we normally used 4 or 5 sectors, 
with each additional sector being ten times longer than the 
preceding one. This technique has enabled us to reduce the 
computational time by up to three orders of magnitude, 
depending on the funnel geometry. Although the grid ex- 
tension can in principle be continued indefinitely, there are 
other factors that limit how far along the jet one can ad- 
vance in practice. Firstly, once the paraboloidal jets become 
highly collimated the required number of grid cells along the 
jet axis increases, and each successive sector becomes more 
expensive than the previous one. Secondly, computational 
errors due to numerical diffusion gradually accumulate in 
the downstream region of the flow and the solution becomes 
progressively less accurate (see Fig. [4]| . 



4 RESULTS 

Models A, B, C and D have geometrical power indices 
a = 1, 3/2, 2 and 3, respectively. Further classification is 
based on the rotation law: models Al-Dl have non-uniform 
rotation, whereas models A2-D2 have uniform rotation. 

Models with different power indices but the same ro- 
tation law show remarkably similar properties. Thus, it is 
sufficient to show only one of them in greater detail. For 
this purpose we selected models CI, C2 and A2. 

Fig. [T] shows the distribution of the Lorentz factor, the 
lab-frame rest-mass density, the poloidal magnetic field, and 
the poloidal electric current for model CI. The top panel 
presents the solution in the innermost grid sector whereas 
the bottom panels show the solution in the second and third 
sectors "glued" together. The density distribution as well 
as the magnetic field lines clearly indicate that the jet de- 
velops a core where the magnetic surfaces become almost 



cylindrical^ The core is produced by the hoop stress of the 
toroidal magnetic field that is wound up in the main body 
of the jet due by the rotation at its base. The ratio of the 
core radius to the jet radius decreases monotonically with 
increasing distance from the source until the core eventually 
becomes unresolved on the grid. The solution then develops 
what looks like an axial line current. (Such behaviour is also 
observed in models A2-D2; see Fig. 1 1 1 p . 

Near the jet boundary, where the angular velocity of 
the magnetic field lines in model CI vanishes, the azimuthal 
magnetic field component is weak and the equilibrium is 
supported in part by the poloidal field. If the flow were 
self-similar, one would have B v oc rj 2 , oc r" 1 , and the 
pressure of the azimuthal field in the main body of the jet 
would eventually become much larger than the pressure of 
the poloidal field in the boundary sheath, leading to a loss of 
force balance. In reality, the flow adjusts through a progres- 
sive compression of the sheath. Consequently a thin layer 
of surface current gradually develops as the distance from 
the source increases. This is why some current lines in the 
bottom-right panel of Fig. [T] appear to terminate at the jet 
surface. 

The most effective acceleration of the CI jet occurs at 
intermediate cylindrical radii, where the angular velocity of 
the field lines reaches a maximum and the poloidal electric 
current flows across the jet (so the (l/c)j p xB,f, force is maxi- 
mized) . Note that this aspect of the jet behaviour could not 
have been studied with the help of self-similar models, in 
which the current lines do not change direction within the 
jet. At comparatively small distances from the inlet bound- 
ary the maximum acceleration occurs at r max ~ 0.5 r^, but 
further downstream r max ~ 0.25 r, . This is explained by a 
more effective collimation in the inner region of the jet than 
at the jet boundary (see discussion following equation 1291 in 
Section [5T|) . 

A careful inspection of the velocity field in the lower- 
right panel of Fig. [T] reveals an additional region of effective 
acceleration near the jet axis for z > 10 3 . This acceleration, 
however, is unphysical as it is caused by numerical diffu- 
sion/dissipation in the core that result from large gradients 
of the flow variables that develop there. The gradual growth 
of errors in this region is clearly seen in Fig. [4] which shows 
the flow constants as functions of $ at various distances from 
the source. Beyond z = 10 4 the errors become unacceptably 
large and this makes further continuation of the solution via 
grid extension meaningless. We note in this connection that, 
even in the absence of exact analytic solutions, the existence 
of flow constants makes the jet problem a very useful one 
for testing RMHD codes and assessing their performance. 

Fig. [5] shows the distribution of the Lorentz factor, the 
lab-frame rest-mass density, the poloidal magnetic field and 
the poloidal electric current for model C2. One can see that 
a core still develops but that a boundary sheath is no longer 
present. This is because the uniform rotation of the magnetic 
field lines in this model ensures an effective generation of 
azimuthal magnetic field all the way up to the jet boundary. 



6 The isodensity contours are, in fact, more strongly collimated 
than the magnetic field lines (or flow s treamlines) , a tr ait already 
identified previously i n nonr elativistic l lShu et al.lll995f) and mod- 
erately relativistic JlH |1996|) MHD disk outflows. 
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Figure 6. Poloidal magnetic field lines (solid) and £ =const coordinate lines (dashed) for models C2 (left panel) and A2 (right panel). 
In all models the magnetic field lines show faster collimation than the coordinate lines. 
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Figure 7. Evolution of the magnetic flux distribution across the jet with distance from the inlet. Left panel: model CI. Middle panel: 
model C2. In both cases the solid line corresponds to tj = 10, the dashed line to r\ = 10 2 , the dash-dotted line to r\ = 10 3 , the dotted 
line to n = 10 4 and the dash-triply-dottcd line to r\ = 10 s . Right panel: model A2. The solid line corresponds to 77 = 1, the dashed line 
to i) = 30, the dash-dotted line to n = 3 X 10 2 , the dotted line to r; = 3 X 10 3 and the dash-triply-dotted line to n = 3 X 10 4 . Note that 
in the conical case we use the spherical coordinate 9 = arctan £ (in radians) rather than the £ coordinate. 
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Figure 8. Fcr (solid line), /u (dashed line) and V (dash-dotted line) along a magnetic field line as a function of cylindrical radius for 
models CI (left panel), C2 (middle panel) and A2 (right panel). 
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Fig. [5] shows the development of an axial line current in this 
solution, a result of the gradual decrease of the core radius 
relative to the jet radius (similar to what is seen in model 
CI). Note, however, that the light cylinder is unresolved at 
the distance where the line current is observed. Thus, what 
looks like a line current could be a smoothly distributed 
current inside the light cylinder. 

Inside the jet the electric current flows inward every- 
where and current closure is achieved via a surface current. 
The radial component of the current peaks near the bound- 
ary, resulting in a higher (l/c)j p x force and a more 
effective plasma acceleration in this region. 

As in the CI solution, the numerical errors in model C2 
grow most rapidly near the jet axis (see Fig. [U, although 
they are somewhat smaller in this case. Moreover, the most 
interesting region of the flow, where the acceleration is most 
effective, is now far from the axis and does not suffer from 
these errors as much as in model CI. This feature is char- 
acteristic not only of models C but of all the other models 
as well. For this reason we decided to focus our attention on 
the models with uniform rotation, A2-D2, and in the rest 
of this section we present results mainly for these solutions. 
This choice is further motivated by the fact that models with 
a non-uniform rotation do not seem to exhibit any signifi- 
cant differences with respect to the uniform-rotation models 
besides those that we have already described. 

Given the results of previous analytical and numerical 
studies, which suggested poor self-collimation of relativistic 
magnetized flows (see references in Section [5.1[) . one could 
have expected the magnetic flux surfaces to almost mirror 
the imposed shape of the jet boundary. However, our re- 
sults indicate that the outflows collimate significantly faster, 
and that this property is manifested not only by jets with 
paraboloidal boundaries but also by the ones that are con- 
fined by a conical wall (see Fig. [3j . Fig. [6] shows the mag- 
netic flux surfaces and the coordinate surfaces £ = const for 
models A2 and C2. In both cases the magnetic flux surfaces 
clearly do not diverge as fast as the coordinate surfaces. 
This effect is further demonstrated by Fig. [7] which shows 
the evolution of the magnetic flux distribution across these 
jets (as well as the jet of model CI) with distance from the 
origin. It is seen that the magnetic flux becomes progres- 
sively more concentrated toward the symmetry axis as the 
flow moves further downstream. 

The left and middle panels of Fig. [5] show the evolution 
of [i, Fa and T along selected magnetic surfaces for models 
CI and C2. For model CI this flux surface is in the mid- 
dle part of the jet, where the flow accelerates most rapidly; 
it encloses ~ 1/3 of the total magnetic flux in the jet. For 
model C2 this surface is near the jet boundary, enclosing 
~ 5/6 of the total magnetic flux in the jet. One can see 
that n remains very nearly constant on the surfaces, indi- 
cating that the flow has reached a steady state and that 
the computational errors that we have described above are 
fairly small. The Lorentz factor at first grows linearly with 
cylindrical radius but then enters an extended domain of log- 
arithmic growth. The linear behaviour was previously found 
in the ma gnetically dominated reg ime of self-similar solu- 
tions (e.g. IVlahakis fc Konig]||2003al ). whereas the logarith- 
mic behaviour was shown to characterize th e acceleration in 
the a symptotic matter-dominated zone (e.g. Begchnan fc Lil 
1 19941 ). The range of Lorentz factors in the solutions derived 



in this paper is evidently too narrow to allow us to probe 
the linear growth regime, but we expect that this could be 
done in our forthcoming paper where we consider higher-r^ 
flows. 

The magnetization function a eventually becomes less 
than 1, signaling a transition to the matter-dominated 
regime. The right panel of Fig. [8] shows the evolution of 
[l, Ta and V along the magnetic flux surface of model A2 
that again encloses ~ 5/6 of the total magnetic flux in the 
jet. This conical jet also exhibits a very effective initial ac- 
celeration and a transition to a matter-dominated regime. In 
this case the growth of the Lorentz factor saturates when it 
reaches F ~ 10, a value that corresponds to an acceleration 
efficiency T//J, of 77%. Although the setup of our conical 
jet model most closely evokes the conical flow geometries 
that have in previous works produced very inefficient ac- 
celerations (see Section [5.1[) . the results displayed in Fig. [8] 
demonstrate that this case is not inherently different from 
the other ones. We discuss the reasons for this in Section [5. II 
It is also worth noting that, while our choice of flux surfaces 
in Fig.[5]was arbitrary, the behaviour along these surfaces is 
fairly representative. To demonstrate this we show in Fig. [9] 
the variation of the same parameters across the jets at large 
distances from the inlet. On can see that the jets are matter- 
dominated throughout the entier cross-section. 

Fig. [10] compares the growth rates of V in models A2- 
D2. The numerical errors in these models are less restric- 
tive than in models Al-Dl and make it possible to extend 
the simulations to larger spatial scales. Each of the plot- 
ted curves corresponds to the magnetic surface near the jet 
boundary that encloses ~ 5/6 of the total magnetic flux. 
The left panel shows V as a function of the cylindrical ra- 
dius normalized by the light- cylinder radius. The most in- 
teresting feature of this figure is the very similar growth of 
r for all models. In fact, up to r ~ 10 — 50 n c the curves for 
models A2, C2 and D2 are almost identical. In model B2 the 
Lorentz factor increases somewhat more slowly. Further in- 
spection reveals another anomaly of model B2 - in contrast 
to the C2 and D2 cases, where the highest Lorentz factor is 
found at the jet boundary, the fastest acceleration in the B2 
solution occurs somewhat off the boundary. The reason for 
these anomalies is not clear but it may have something to do 
with the curvature of magnetic field lines - given the lower 
value of the power-law index a, model B2 retains a higher 
curvature at larger radii than models C2 and D2. The rea- 
son why the model D curve is significantly shorter than the 
other is that the strong collimation of the jet rapidly renders 
the computation prohibitively expensive in this case. 

Since more rapidly collimated jets reach the same cylin- 
drical radius at a larger distance from the source, the similar 
growth rates of the Lorentz factor with cylindrical radius im- 
ply a faster growth with spherical radius for less collimated 
jets. This is exactly what we see in the right panel of Fig. [10] 
- the conical jet of the A2 model reaches a Lorentz factor 
of 10 at a distance from the origin that is almost 100 times 
shorter than that of the paraboloidal jet of model C2. 

Fig. [TTJ compares the magnitudes of the different mag- 
netic field components in models A2-C2 near the far end of 
the jet (r) = 10 3 ). At this distance the jet radius is almost 10 3 
larger than the light-cylinder radius and one would expect 
the azimuthal component of the magnetic field to dominate. 
Indeed, this is what is observed. On these scales the light 
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Figure 9. Distribution of fj, (solid line), T (dashed line) and Tcr (dash-dotted line) across the jet for models CI (left panel) and C2 
(middle panel) at r\ = 10 5 , and for model A2 (right panel) at r\ = 2 X 10 3 . 
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Figure 10. Lorentz factor along the jet boundary as a function of the cylindrical radius r (left panel) and the spherical radius R (right 
panel): model A2 - solid line; model B2 - dashed line; model C2 — dash-dotted line; model D2 — dotted line. Note that, although ro and 
r; c differ from model to model, Rg = 1 in all cases. 



cylinder is no longer resolved on the computational grid, 
which explains why the azimuthal field component exceeds 
the poloidal components even near the symmetry axis. The 
fact that does not vanish shows that the solution devel- 
ops a core of high electric current density, and this core is 
also unresolved in our solution. The development of an ax- 
ial line current in model C2 is shown in Fig. [5] very similar 
results are found also for the other models. 

In an attempt to clarify the role of the imposed wall 
on the acceleration and collimation of the simulated jets, we 
performed additional simulations of model A2, correspond- 
ing to higher opening half-angles of the wall: C = 0.4, 0.8, 
and 7r/2, where the last case represents an unconfined out- 
flow[H The results are presented in Figs. \T%\ and 1131 which 



Our original setup for this model corresponded to 8 C = 0.2. 
We keep A9 (the cell size in the polar direction) the same for all 
models by increasing the number of cells for larger 8 C . 



should be compared with the corresponding plots in the right 
panels of Figs. [THS] Not surprisingly, we find that the effi- 
ciency of magnetic acceleration is lower for weaker external 
collimation. This is particularly apparent in the unconfined 
case (0 C = 7r/2), where the flow remains Poynting-dominated 
in the equatorial sector at all distances. Nevertheless, even 
in the equatorial region of this extreme case the acceleration 
efficiency Voo/fi is ~ 40%. Furthermore, in the polar zone 
the efficiency is high, similarly to the small-# c case, with 
all models exhibiting a transition to the matter-dominated 
regime. In this sector self-collimation is strong as well (see 
Fig. 1 1 3p . The behaviour of ou r unco nfined jet model is con- 
sistent with the results of Q (|l996t ). who showed (for both 
the nonrelativistic and modestly relativistic regimes) that 
a collimated axial jet can form from an initially spherical 
MHD wind. Those results were, however, derived on the as- 
sumption that the magnetic field is purely azimuthal and 
the velocity field purely poloidal from the start. The good 
collimation exhibited by all our jet models in the vicinity of 
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Figure 11. Variation of the magnetic field components across the jet at jj = 6 X 10 3 in models CI (left panel), C2 (middle panel) and 
A2 (right panel): log 10 \B^\ - solid line, log 10 B* 1 - dashed line, log 10 \B^\ - dash-dotted line. 




0.1 0.2 0.3 0.4 0.2 0.4 0.6 0.8 0.5 1 1.5 

see 

Figure 12. Distribution of fj, (solid line), T (dashed line) and IV (dash-dotted line) across the jet for the cases of a conical boundary 
with an opening half-angle 9 C = 0.4 (left panel), 0.8 (middle panel), and n/2 (right panel, an unconfincd flow) at r\ = 2 X 10 3 . These are 
to be compared with the results for 6 C = 0.2 shown in the right panel of Fig. [9] 
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Figure 13. Left panel: Variation of To (solid line), fi (dashed line) and T (dash-dotted line) along a magnetic field line for the case of a 
conical boundary with 8 C = it/2 (i.e. an unconfincd flow). This is to be compared with the results for 8 C = 0.2 shown in the right panel 
of Fig. [8] Right panel: Evolution of the magnetic flux distribution across the unconfined jet. The solid line corresponds to Tj = 1, the 
dashed line to t) = 30, the dash-dotted line to r\ = 3 X 10 2 , the dotted line to r\ = 3 X 10 3 and the dash-triply-dotted line to n = 3 X 10 4 . 
This is to be compared with the results for 9 C = 0.2 shown in the right panel of Fig. [7] 
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the axis can be attributed to the action of magnetic hoop 
stresses associated with a poloidal current flow, as discussed 
in §0 



5 DISCUSSION 

5.1 Theoretical aspects of the problem 

Over the years there have been persistent doubts in the lit- 
erature regarding the ability of magnetic forces to accelerate 
flows to relativistic speeds. In particular, several published 
studies have concluded that MHD acceleration of relativis- 
tic flows is inherently inefficient. This conclusion, however, 
is erroneous and can be attributed to the adoption of a 
conical (split-monopole) flow geometry in these studies. For 
example, in the work of iMichell (|l969h a simplified conical 
geometry was used in which the full system of relativistic 
MHD equati o ns w as not satisfied, whereas the results of 
iBeskin et ail (ll998T l were based on a perturbative analysis 
around a quasi-conical flow. The conical flow geometry is 
unfavorable for acceleration for the following reason. Well 
outside the light cylinder, where rQ 3> and v ~ c, equa- 
tions (J23J) and ([24} imply 



are not conical but rather paraboloidal (see Figs. [6] and [7}, 
so B p r 2 is not constant. Fig. 1141 shows the evolution of the 
function 



rB^ = --QB p r 2 . 

c 



From this equation and equation (|16p one finds that 



I = -\o.B v r 2 . 
2 



(27) 



(28) 



where B p is the magnitude of the poloidal magnetic field. 
If the magnetic surfaces are conical then B p oc r~ 2 , and 
thus the poloidal electric current flows parallel to the mag- 
netic field lines. In this case the component of the Lorentz 
force along the poloidal magnetic field lines, (l/c)j p xB^, 
simply vanishes. More general treatments of the problem, 
based on exact semi-analytic solutions for axisymmetric, 
highly magnetized, stea dy outflows under the assumption 
of rad ial self-s imilarity (|Li et al.l 1 19921 ; IVlahakis fc KonigH 
l2003al lbl. I2004T ). have demonstrated that magnetic acceler- 
ation in non-conical geometries can be quite efficient, typi- 
cally resulting in a rough asymptotic equipartition between 
the Poynting and matter energy fluxes. A similar conclusion 
was reached o n the basis of a perturbati ve analysis around a 
parabolic flow (Bes kin fc Nokhrina 2006). These results have 
indicated that the correct paradigm should, in fact, be that 
magnetic acceleration is generally a rather efficient mecha- 
nism for producing relativistic flows. 

In this paper we have for the first time verified the 
proposed paradigm by means of numerical simulations of 
highly magnetized, relativistic flows. We have focused on 
the parameter regime that is most relevant to AGN jets. In 
a future paper (Komissarov et al., in preparation) we will 
present additional simulations that will demonstrate that 
this paradigm also applies to flows with terminal Lorentz 
factors that are as high as those inferred in gamma-ray burst 
sources. 

One of the interesting outcomes of this study is the 
highly effective acceleration even in the case where the shape 
of the outer boundary is conical. Although the acceleration 
efficiency in conical steady flows is - as explained above - 
tiny, our results show that magnetic surfaces of conical jets 
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along a typical magnetic surface for models C2 and A2. It 
is seen that in both cases S undergoes a significant decrease 
with distance from the source. In fact, this decrease is faster 
in the conical-boundary model, which is reflected in the more 
rapid acceleration in this case (see Fig. II Of) . The direct re- 
lationship between the function S and the acceleration ef- 
ficiency can be readily shown by combining equations (|17[1 
and (1281) to obtain 



To- = (*fi 2 / 47I " 2fcc3 )5 <xS. 

Since S depends on the shape of the flow, the latter 
relation brings out the importance of the trans-field force 
balance and the connection between acceleration and col- 
limation. If the poloidal magnetic field is almost uniformly 
distributed across the jet then S ~ 1; this is the case near the 
inlet boundary. However, due to the collimation, the poloidal 
magnetic flux becomes concentrated near the rotation axis, 
forming a cylindrical core and causing S to decrease with 
increasing r (see Fig. I15[) . 

Our jet solutions are characterized by a high magnetic- 
to-kinetic energy conversion efficiency, but in the final states 
that we obtain the kinetic and Poynting energy fluxes are 
still of comparable magnit ude, as previously found i n the 
self-similar solutions (e.g. Vlahakis & Konigl 2003a). We 
note in this connection that iBegelman fc Lil (|l994T ) conjec- 
tured on the basis of their asymptotic analysis that relativis- 
tic MHD outflows with a finite total magnetic flux will tend 
to convert their entire magnetic energy into kinetic energy at 
large distances from the origin. Although our solutions are 
terminated at a finite distance and thus we cannot exclude 
the possibility that eventually nearly all the magnetic energy 
will be extracted, it appears that this is unlikely to occur 
on astrophysically realistic scales. Further more, our discus- 
sion a bove suggests that the reason given by Begclman fc Lil 
|1994| ) for the limited conversion efficiency of self-similar 
flows, namely, that it may be hard to meaningfully decrease 
the quantity nr 2 B p in a medium that contains an infinite 
magnetic flux it, is not correct, since the relevant quantity 
is 5 = -Kr 2 B p /^i (equation I29p . which remains finite even in 
the self-similar case. 

Our results show efficient self-collimation, contrary to 
some claims in the literature and certain seemingly simi- 
lar investigations that reached the opposite conclusion. As 
the term "self collimation" is sometimes misunderstood and 
confused with "self confinement", we start by clarifying its 
meaning. In general, a hydromagnetic jet cannot be "self 
confined" for the simple reason that the "outward" effect of 
magnetic pressure always overcomes the "inward" effect of 
magnetic ten sion in MHD equilibria (see e.g. Section 5.3 in 
|Parker|[l979h . There must always be an external medium 
that confines the system at its boundaries. (In our case 
the geometry of the wall translates into an effective am- 
bient pressure distribution for any given flow problem; see 
Fig. 1171 below.) It is, however, possible for a jet that carries 
an axial current to have the magnetic hoop stress collimate 
the innermost (current-carrying) streamlines relative to the 
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outer regions of the flow. This is what is generally meant by 
"self collimation." The outer streamlines may or may not 
be collimated, depending on the glob al current flow and th e 
confining pressure distribution (e.g. iBegelman fc Lilll994h . 
The behaviour of the innermost streamlines is determined 
mostly by the axial current distribution. By showing that 
these streamlines collimate faster than the imposed bound- 
ary, we were able to demonstrate the action of the "self col- 
limation" process in our jet models, including the limiting 
case of an unconfined flow (Fig. I13|l . The different cases that 
we examined served to illustrate how the transition from the 
inner region (small cylindrical radius R) to the outer region 
(large R) occurs for different effective confining pressure dis- 
tributions. 

Self-collimation in the super-Alfvenic regime of magne- 
tized outflows is the result of the (l/c)j p xB^ force in the 
trans-field direction. In relativistic flows, the effect of this 
force is almost completely countered by the electric force, re- 
sulting in slower collimation compared to the nonrelativistic 
case (where the electric force is negligible). The asymptotic 
form of the trans-field equation in the highly relativistic limit 
is 
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Figure 14. Evolution of the function S = rrB v r 2 along the 
magnetic surfaces of models A2 (solid line) and C2 (dashed line). 
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(see equation 16 in Vlahakis 2004). From this equation it 
follows that the radius o f curvature of poloidal field li nes 
is TZ ~ T 2 r. This fact led iTsinganos fc Bogovalovl i|2002T ) to 
propose a two-component outflow model (central jet and sur- 
rounding disk wind) as a way of explaining the collimation 
of relativistic jets. Howev er, as the self-similar solutions of 
IVlahakis fc Konigll \2004 ) as well as the present simulations 
show, self-collimation is still possible. In fact, it remains pos- 
sible in flows with even h igher asymptotic Lorentz factors 
I Vlahakis fc Konigll [2003al ') . Although for V > 1 the colli- 
mation is indeed slow, it is more efficient near the source, 
where the flow is not yet highly relativistic. 

iBogovalovl (|200lh solved a similar problem (using time- 
dependent equations near the central source and steady- 
state equations further out). Although the setup in that pa- 
per is similar to our case A2, the conclusions are different (in- 
efficient collimation and therefore less efficient acceleration 
compared to our solution). There is, however, an important 
difference in the two setups. In Bogovalov's (2001) paper 
the poloidal velocity at the inlet boundary is v po ~ 0.87 c, 
corresponding to a Lorentz factor (including the azimuthal 
velocity) significantly higher than in our simulations. As ex- 
plained above, a high Lorentz factor leads to a large TZ/r 
(~ r 2 ). Another difference between the two works is that 
we are able to follow the flow to larger distances and hence 
to a smaller-cr regime: in fact, in some cases (A2 and B2) 
our solutions extend all the way to the asymptotic regime, 
where the acceleration ceases and the Lorentz factor satu- 
rates to a constant value. Related to the above discussion is 
the fact that the mass and magnetic flux in our jet solutions 
are not "uncomfortably low" (as th ey were acknowledged to 
be in Bogovalov's 2001 solution; see ITsinganos fc Bogovalovl 



2002). In our solutions all the magnetic flux and all the out- 
flowing mass are effectively collimated[f] 

The preceding discussion of the magnitude of the radius 
of curvature implicitly assumes that TZ is positive. Accord- 
ing to equation (|30|) . the sign of TZ depends on the sign 
of the three quantities V|J| • V*, Vr • V* and Vr • V*. 
The term Vr ■ V*& always corresponds to decollimation and 
is important onl y in the matte r-dominated flow regime far 
from the source ( Vlahakis 2004). We can ignore this term in 
the main acceleration region where the flow is still magneti- 
cally dominated. The j p B p < current-carrying regime (in 
which V|I| ■ > 0) contributes to positive TZ and thus to 
collimation, whereas the return-current regime j p ■ B p > 
promotes decollimation ((Okamot o 2003). However, the sign 
of TZ also depends on the gradient of F, a manifestation of 
the electric force. It is possible to have TZ > even in the 
return-current regime provided that the Lorentz factor de- 
creases with increasing cylindrical radius (Vr ■ V*& < 0). In 
this case the electric force, which is directed toward the axis, 
dominates over the magnetic force, which points away from 
the axis, leading to collimation0 The net effect of the total 
electromagnetic force depends on the gradient of |/|/F, and 
collimation is possible if this quantity increases on moving 
acros s the field lines away from the polar axis (see also |lj 
1996). In agreement with this analysis, the positive value of 
TZ in our models Al-Dl (which contain a current-carrying 
region near the axis and a return-current region near the 
outer boundary) requires the Lorentz factor to decrease with 
cylindrical radius (see e.g. Fig. [1} . The decrease in F as the 
outer wall is approached is consistent (by equation I17|l with 



8 The mass-loss rate between the axis and a particular field line 
\I/ =const is M = 2 J" k(W)<M!. Since fc(\I/) is practically constant 
(see Fig.3{, M oc f and the distribution of mass- loss rate across 
the jet can be deduced from the behaviour of \P(r) in Fig. 1151 

9 This behaviour was manif ested also by the return-cur rent self- 
similar solution presented in I Vlahakis fc K onigl (20 03al l. 
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the reduction in the electromagnetic acceleration brought 
about by the imposition of the boundary condition = 
at £ = £j (see equation I21[) . 

To summarize the discussion on the collimation effi- 
ciency of relativistic jets, we have argued that collimation 
is possible in accelerating flows where the Lorentz factor 
ranges from To ~ 1 near the source to a high asymptotic 
(subscript oo) value r x . Our choice of v P0 (= 0.5 c) at the 
inlet boundary and of the initial value of v v for the funnel 
plasma (= 0.7 c) allows the flow to relax to a collimated 
steady state with a high Toe. The sign of the curvature ra- 
dius is positive even in the return-current regime because the 
Lorentz factor decreases sufficiently rapidly with r across the 
jet. 

According to the asymptotic analysis of 



iHevvaerts fc Norman (1 1989h and its relativistic gener- 



alization by Chiueh et alj (Il99ll ). the formation of a 

way to have a non-singular 
3 (* = 0) = 0. Although our 



cylindrical core is the only 
current near the axis, i.e. I 
numerical results do not include the far-asymptotic state, 
the "solvability condition" J r oc(5 , )/r oo ( , I') =const (a direct 
consequence of equation [30] in the limit r/TZ — > + ) is 
roughly satisfied; see Fig. If 61 

The profile of fi(£) on the inlet boundary is used in our 
simulations as a proxy for the current distribution (see equa- 
tion [28]). However, in the cases Af-Dl, where fi vanishes 
on the outer boundary, I nevertheless remains finite there. 
(Note that equation l28l holds only for rQ ^> v^, which is not 
satisfied when Q vanishes.) This is a reflection of a general 
property of ideal MHD flows: in the super- Alfvenic regime 
the azimuthal magnetic field cannot vanisfF°l and the cur- 
rent lines close in a current sheet. In this respect our choice 
of an outer "wall" captures a basic physical aspect of a real 
boundary between a jet and an unmagnetized environment. 

Since we have constructed our numerical jet models us- 
ing time-dependent simulations, we have partly addressed 
the issue of stability to axisymmetric perturbations. Our 



10 As explained in IVlahakisI ||2004| ') . the invariance of B 2 — 
E 2 yields B^/B 2 > (r 2 /rf c ) — 1. In highly magnetized flows 
the Alfven surface almost coincides with the light surface, so 
(r 2 /rf c ) - 1 > in the super- Alfvenic regime. 



results indicate that, within our model setup, the jets are 
stable. Although the imposed solid outer boundary could 
have some stabilizing effect even on these modes, the ex- 
tent of this influence is unclear, although in any case we 
do not expect the outer wall alone to prevent the growth 
of the pinch mode in the jet interior. Linear stability anal- 
yses of "Z-pinch" equilibria (magnetostatic cylindrical con- 
figurations with a weak uniform poloidal field and a strong 
toroidal field) have been used to argue that the kink mode 
would destroy the co ncentric configurations of astrophysical 
jets (jBegel man 1998). While our simulations cannot address 
the possible effect of nonaxisymmetric modes, we neverthe- 
less note that, in contrast to "Z-pinch" configurations, our 
jets show strong transverse gradients of poloidal velocity 
and magnetic field, w hich may be stabilizing factors (e.g. 
I Anderson et al.ll2006bl ). 



5.2 Application to AGN Jets 

The initial energy-to-mass flux ratio of jets in our simula- 
tions yields an upper limit on the terminal Lorentz factor 
Too = /i < 16. This is consistent with the mean values in- 
ferred in AGN jets (see Section [TJ. In order to make further 
comparisons of our numerical models with observations we 
need to select suitable dimensional scales. The key scale in 
the problem of magnetic acceleration is the light cylinder (or 
the Alfven surface) radius, r; c . If the jets are launched by a 
rapidly rotating black hole in the center of an AGN then 

r ic ~ 4r 9 = 6 x 10 13 (M/10 8 Mq) cm, 

where r g = GM/c 2 . In this estimate we assume that the 
angular velocity of the magnetic field lines is half of that of a 
maximally rotating (rotation parameter a ~ 1 .0) black hole. 
According to the results shown in Fig. [5] the jets enter the 
matter-dominated regime, where an equipartition between 
Poynting and matter energy fluxes (a ~ 1) is established, at 
a cylindrical radius 

r eq ~ 30r ic ~ 2 x f0 15 (Af/1 8 M ) cm, 
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Figure 16. Distribution of I/T across the computed jets. Left panel: Model CI at r\ = 80 (solid line), r\ = 800 (dashed line) and n = 8000 
(dash-dotted line); Middle panel: Model C2 at rj = 80 (solid line), n = 800 (dashed line) and r\ = 8000 (dash-dotted line); Right panel: 
Model A2 at r\ = 20 (solid line), n = 200 (dashed line) and n = 20000 (dash-dotted line). 



more or less independently of the details of jet collimation. 
The corresponding distance from the black hole is 




where 0j is the jet opening half-angle. If blazar flux vari- 
ability is associated with the propagation of strong shocks 
within the jet then we can expect this behaviour to origi- 
nate on scales >R eq . When our simulated jets reach R ~ 
10 Req, their characteristic Lorentz factor becomes ~ 10. 
These properties of the extended magnetic acceleration re- 
gion are in very good agreement with the observational in- 
ferences summarized in Section [T] (in particular, the lack 
of bulk-Comptoni zation spectral features in blazar jets; see 
ISikora et al]|2005h . 

Taking the characteristic initial radius poloidal mag- 
netic field of a black hole-launched jet to be ro = r g and 
Bo = 10 5 G, respectively, the mass-loss rate and total lumi- 
nosity of the model jets scale as 




respectively, where k, \P and \x are the mean values of the 
dimensionless flow constants shown in Fig. [4] 

From the theory of black-hole magnetospheres (e.g. 
iBlandford fc Znaiekl fl977l ) it follows that, at their base, 
black-hole jets are highly magnetically dominated, so that 
the energy per particle greatly exceeds the Lorentz factor 
inferred from observations of AGN jets. This difficulty can 
be overcome if there are other ways of injecting particles 
into AGN jets in addition to pair cascades. It is conceiv- 
able that a sufficient supply of particles is provided by the 
winds of stars that lie in the paths of the jets as they make 
their way out of the galact ic nuclei (e.g. iKomissarovl 1 1994 ; 
Irlubbard fc Blackman|[2006l ) . In fact, the injection rate could 
be high enough to explain the observed deceleration of weak 
FR-I jets down to subrelativistic speeds. In addition some 
mass is unavoidably entrained through the interface between 



the jet and the confining external medium via diffusive pro- 
cesses. Another possibility is that the speed of AGN jets is 
controlled by the dissipation at oblique shocks developing as 
a result of an interaction with an unsteady external flow or 
due to various instabilities (|McKinnev |2006). Compton-drag 
interaction with the ambient rad iation field may also be a 
factor (e.g. lMelia fc Konigl|[l989l ). 

The problem of low initial mass loading might be cir- 
cumvented if the bulk of the relativistic outflow compo- 
nent in_Ja£t_OTj|£iiiate^ (see, 
e.g. iGhosh fc Abramowic j 1 19971 ; iLivio et all Il999h . which 
is the scenario adopted in the semi-analytic models of 
IVlahakis fc Konigll (|2004l ). Although the chosen distribu- 
tions of angular velocity, magnetic field and density at the 
inlet boundary do not formally match a Keplerian disk, our 
solutions can be interpreted in the context of disk-driven 
outflows. Taking £!k(j"o) to be the Keplerian angular veloc- 
ity at the reference distance ro, we find that 



( M ) 




^3/2 


f2 K (r ) 











and that a ~ 1 is reached at a distance 

1017 fa^) (Iff 



The mass-loss rate and jet power are 




respectively. The striking differences in the velocity profile 
between simulated jets with solid-body rotation and with 
differential rotation (see Figs. [T] and [2] for models CI and 
C2, respectively) suggests (taking into account relativistic 
beaming effects) that it might be possible to distinguish be- 
tween jets launched directly from a black hole and those 
that emanate from the surface of an accretion disk in cases 
where the transverse structure of the jet can be resolved. 
One would, however, need to verify that these differences 
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Figure 17. Ambient pressure distribution required for the col- 
limation of the computed jets. Model A2 - solid line; model B2 
- dashed line; model C2 - dash-dotted line; model D2 - dotted 
line. 



remain noticeable for more realistic disk rotation laws and 
surface field distributions. 

The funnel shapes in our simulations were chosen 
merely for the purpose of studying the effects of overall 
flow collimation on its magnetic acceleration. However, from 
our steady-state solutions we can infer the effective external 
force (normal to the jet surface) that is required to provide 
the collimation imposed by our choice of the outer-boundary 
shape. AGN jets could be confined by a variety of forces, in- 
cluding, for example, the thermal pressure of an ambient 
gas distribution, the ram pressure of a wind from the outer 
regions of the nuclear disk, and the stress of a magnetic 
field anchored in the disk (and possibly embedded in a disk 
outflow). Fig. ll7l shows the effective pressure deduced in this 
way (pext = b 2 /8ir) as a function of spherical radius for mod- 
els A2-D2. Although none of the curves is an exact power 
law of the form jw oc R~ a , it is nevertheless informative 
to calculate mean power-law indices. We find a « 3.5, 2, 1.6 
and 1.1 for models A2, B2, C2 and D2, respectively. The 
models are thus seen to cover a wide range of behaviours. 
As expected, the more highly collimated funnel geometries 
correspond to the less steeply declining effective pressure 
distributions. The largest indices might correspond to con- 
finement by a wind. For example, in a spherical wind of poly- 
tropic index 5/3, the thermal pressure scales as R~ 10 ^ 3 and 
the ram pressure as R~ 5 / 2 . Thus, a disk wind that assumes 
a nearly spherical geometry as it propagates away from the 
disk surface could effectively confine a relativistic jet with 
a nearly conical outer boundar y. We note in this connec- 
tion that Bcgc lman fc Lil l| 19941 ) suggested that collimated 
jets necessarily correspond to confining pressure distribu- 
tions with a<2. Our results indicate that higher values of 
the effective power-law index are also possible. 

Although we have focused in this paper on AGN jets, it 
is interesting to note that relativistic outflows with terminal 
Lorentz factors as high as ~ 10 have also been inferred in 



Galactic X-ray binary source s, which comprise both black 
holes and neutron stars (e.g. iFender et alj 12004). and that 
arguments have been advanced in support of the possibility 
that the mean Lorentz factors in these sources are compara - 
ble to those estimated in AGNs (jMiller- Jones et al.ll2006)l) . 
The magnetic acceleration mechanism discussed in this pa- 
per is als o a likely candida te for the driving of X-ray binary 
jets (e.g. iLivio et al.ll2003l ). However, even if these jets are 
similar to those in AGNs, as of now the latter remain the 
best targets for observations that could test and constrain 
the model. 



6 CONCLUSION 

This paper presents the results of special-relativistic, ideal- 
MHD numerical simulations of AGN jets. The numerical 
code employed in these simulations was specifically designed 
for this task. In contrast to most previous numerical schemes 
that modeled relativistic MHD jets, our code does not re- 
quire a large artificial viscosity for numerical stability and 
is well suited for studies of two-dimensional stationary flows 
that are aligned with the computational grid. To avoid nu- 
merical dissipation induced at the interface with an ambient 
medium, we have simplified the calculation by taking the 
flow to be confined by a solid wall. We took the shape of 
the wall to be either a paraboloid of revolution or a cone 
and used corresponding elliptical or spherical coordinates to 
optimize the resolution as well as the alignment of the flow 
with the computational grid. In addition, we implemented 
a grid-extension method that allowed us to follow the flow 
out to scales ~ 10 4 — 10 6 times that of the inlet boundary, 
which was crucial to our ability to study the inherently ex- 
tended nature of MHD acceleration to high Lorentz factors. 
To ensure the self-consistency of this procedure, we derived 
the condition for the Mach cone of the fast-magnetosonic 
waves to point outward at the boundary between a given 
grid sector and the successive one, and we verified that this 
condition is satisfied at each of the relevant sector interfaces. 

Our carefully designed numerical scheme has enabled us 
to simulate, for the first time, the magnetic acceleration and 
collimation of relativistic jets to completion. In particular, 
we have found that initially Poynting flux-dominated jets 
can be effectively accelerated to high bulk Lorentz factors 
with an efficiency (defined as the ratio of the final kinetic 
energy flux to the total energy flux) > 50%. As expected 
from previous semi-analytic (radially self-similar) solutions 
for steady-state flows, the acceleration process is spatially 
extended. We have found that our simulated jets invariably 
settle to a steady state, which suggests (although we did 
not explore this issue explicitly) that the resulting flow con- 
figurations are not inherently unstable (at least not to ax- 
isymmetric perturbations - although it is conceivable that 
the imposed rigid wall has a stabilizing influence in this re- 
gard - and excluding by design any effects of a direct in- 
teraction with an ambient medium). The properties of the 
derived final configurations were found to be qualitatively 
ve ry similar to those of the self-similar AGN jet solutions 
of IVlahakis fc Konigll l|2004l) and to not depend sensitively 
on either the imposed shape of the outer boundary or on 
the distribution of the injected poloidal current at the inlet 
boundary. (We explored boundaries with scalings, in cylin- 
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drical coordinates, ranging from z <x r to z <x r 3 , and current 
distributions that either closed within the volume of the jet 
or on its outer boundary.) 

We provided a physical explanation of the basic acceler- 
ation process and of the variations in the detailed behaviour 
among the different flow configurations that we simulated. 
We argued that the robustness of the acceleration process 
can be attributed to the fact that the bulk of the outflow 
initially follows paraboloidal trajectories, including the case 
of a conical outer boundary. We highlighted the connec- 
tion between the collimation of the flow, which is mani- 
fested in the curved streamlines, and the acceleration pro- 
cess. The collimation in the current-carrying regime is es- 
sentially due to magnetic hoop stress associated with the 
azimuthal magnetic- field component B^. The collimation 
induces a reduction in the magnitude of r 2 B p (where B p 
is the poloidal field component) along the poloidal stream- 
lines, which corresponds to a decrease in the Poynting flux 
along the flow and therefore results in acceleration (driven 
by the gradient of the magnetic pressure associated with 
B$). Previous claims in the literature that magnetic accel- 
eration of relativistic flows is inefficient were all based on the 
assumption that the streamlines have a split-monopole ge- 
ometry (or very nearly so), which is a singular case in which 
by fiat r 2 B p remains constant (or close to a constant) along 
the flow. 

Our solutions also revealed that the collimation effi- 
ciency of relativistic jets can be high if they are acceler- 
ated from an initial Lorentz factor To ~ 1. We argued that 
published results in the literature that claimed otherwise in 
fact had a significantly higher To. Once the flow attains a 
high Lorentz factor the collimation process slows down on 
account of the increased inertia and of the growth of the 
electric force, which almost completely balances the trans- 
verse magnetic force. Nevertheless, the current-carrying cen- 
tral region of our simulated jets collimates much faster than 
the imposed boundaries and attains a cylindrical shape by 
the time the terminal Lorentz factor is attained, again in 
full agreement with the semi-analytic self-similar solutions 
(which also assumed To ~ 1). In simulated outflows where 
the current returns through the jet we found that the flow 
is effectively collimated also in the outer, return-current re- 
gion (in this case by the electric force, which dominates the 
transverse magnetic force that acts to decollimate the flow 
in this regime). The efficient electromagnetic collimation in 
all of our computed jet models is evidently the reason why 
the presence of a rigid outer boundary does not induce rec- 
ollimation shocks in the outflow even for the most rapidly 
converging wall shape (r a z 3 ). 

In validating the basic features of the simplified semi- 
analytic solutions, our numerical results go a long way to- 
ward establishing an "MHD acceleration and collimation 
paradigm" for relativistic astrophysical jets. In this con- 
tribution we applied this model to AGN jets, for which 
there is already significant observational evidence of ex- 
tended, > 0.1 pc-scale acceleration (possibly continuing to 
~ 1— 10 pc) to Lorentz factors > 10. We demonstrated that, 
for plausible physical parameters, our simulated jets can 
repro duce these observations (see also IVlahakis fc Konigll 
120041 '). We noted that these results could potentially apply 
also to the jets observed in Galactic X-ray binary sources. 
In a future publication we will present simulation results 



for even higher initial magnetizations that could be used to 
model the ultrarelativistic jets in gamma-ray burst sources. 

Future investigations will need to examine some of the 
complicating factors left out of the present study. In partic- 
ular, the simplifying wall boundary condition will have to 
be replaced by the interaction of the outflow with its en- 
vironment. Furthermore, the questions of stability and en- 
ergy dissipation will have to be addressed. In fact, given the 
potentially destructive role of the kink mode, fully 3D rel- 
ativistic simulations will ultimately need to be carried out. 
In addition, the response of the jet to nonsteady conditions 
at the source - a potential cause of observed variability - 
will be an interesting topic for study. 
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APPENDIX A: ELLIPTIC COORDINATES 

We assume that the jet boundary satisfies the power law 

z oc r a , (Al) 

where {r, z} are cylindrical coordinates. Then the condition 
that the jet boundary be a coordinate surface suggests we 
choose 



£ = rz 



-l/a 



(A2) 



as one of the spatial coordinates. The other coordinate, r), is 
defined in such a way that the coordinate system becomes 
orthogonal. The orthogonality condition 



V£ ■ Vr/ = 

leads to a PDE for r\ 



drj 1 r drj 
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r a z dz 



0. 



(A3) 



(A4) 



which allows separable solutions. The requirement rj = for 
(r, z) = (0, 0) leads to 



2 r 2 

rj = h z . 

a 



(A5) 



Thus, the rj coordinate lines are ellipses with semi-axes rj 
and y/arj. The remaining spatial coordinate is the usual az- 
imuthal angle <j>. 

Conversion from elliptical to cylindrical coordinates in- 
volves solving a transcendental equation for z: 



z 2 + ^z 2 ' a 
a 



V = . 



(A6) 
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In the general case this equation has no analytic solutions, 
but for certain values of the power-law index a it reduces to 
simpler equations: 



= 3/2, 


3 | 2 2.2 2 

y + g£ v -v = 


0, where y — z 3 ^ 2 ; 


(A7) 


= 1, 


z 2 +£ 2 z 2 - V 2 = 0; 




(A8) 


= 2, 


z 2 + L z -r l 2 = 0; 




(A9) 


= 3, 


? 2 

3 , ? 2 n 

y + -g y-n =o, 


where y = z 2//3 . 


(A10) 



The metric tensor in these coordinates is diagonal with 
components 



fee 



9vv 



D 



2 2 

a rf 
D 



94,4. = r , 
Qtt = -1 , 

where D — a 2 z 2 + r 2 . Its determinant is 



9 



a 4 r 2 v 2 z 2(l+a)/a 

D 2 



(All) 

(A12) 

(A13) 
(A14) 

(A15) 



The non-vanishing derivatives of these components are 
S€5,« = ~j^2a 2 rz (3+2a)/a [2a 2 z 2 + r 2 (l + a)] , (A16) 



1 o if ,n (l + 2a)/a 2 

ffw.S = ^3 2a (a - l)rz k rj . 

^ o 3 / 1 \2 2 2 

ffw.i? = ~Jp 2a (a-l) r z r), 

1 r, 2 (l + 2a)/a 

94>4>, i = -J 2a rz k 

1, 2 

500, »j = -p2ar ??. 



(A17) 
(A18) 
(A19) 
(A20) 
(A21) 



APPENDIX B: FAST MAGNETOSONIC WAVES 

Suppose that we study an axisymmetric magnetosonic dis- 
turbance whose wavevector is 



k — k ( — cos i?fc + 4> x — sin 
Its frequency is u> — kV^k), with V satisfying 



(Bl) 



B 2 c 2 

±jp y^ s 



cos dh 



4nw/c 2 

— (TV — Up cost?*,) 2 



— cos Vk H 75- 

r; c c Bp 



V 2 



B2 - E ' 2 f i-i) + o 2 



Anw/c 2 



(TV — Up cosi9 fe ) 
1 -V 2 /c 2 



]=0 (B2) 



(see Appendix C in lVlahakis fc Koniglll2003ah . The distur- 
bance travels with a group velocity 



k k 

v„ = V k tu = -V + <b x -V' , 



where V' = dV/d-ff^- The group velocity makes an angle # s 
with respect to the poloidal flow velocity, where 



tantf 9 = " 3 = ' a*^' „* . (B4) 

The envelope of the family of such disturbances (whose tra- 
jectories are defined by the angle #&) constitutes the Mach 
cone of fast-magnetosonic waves at any given point in the 
flow. It is given by combining equation ()B4[) and the condi- 
tion 

d ( Psini9 fc + P'costf fc \ _ 
d$k \Vcostfk - V'smfikJ ~ 

After some manipulation, equation IB5I yields V = O ri 
Thus, the fast Mach cone corresponds to a particular com- 
bination ($k , $9) of the angles •dk and i? 9 that satisfies 
tani9 9 = — coti9fe and Vifik) = 0. Using equation (|B2[) . 
these conditions yield 



(B5) 



B - E l 

A-kw/c 2 



+ ■ 



l/c 2 



B 2 p/u 2 v 
Anw/c 2 



(B6) 



The requirement that the Mach cone points outward at 
a surface r\ =const is that the angle arcsin \ fj-v p /v p \ between 
the surface and the poloidal flow velocity exceeds or \rj ■ 
Vp/v p \ 2 > sin 2 {>g. Using equation (|B6[) and B p /v p — B n /v**, 
this last inequality can be transformed into equation 



(B3) 
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